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<^ | Abstract 

\o : A dripping faucet system is simulated. We numerically show that a hanging drop generally 



has many equilibrium shapes but only one among them is stable. By taking a stable equilibrium 
shape as an initial state, a simulation of dynamics is done, for which we present a new algorithm 
■ based on Lagrangian description. The shape of a drop falling from a faucet obtained by the present 

algorithm agrees quite well with experimental observations. Long-term behavior of the simulation 
can reproduce period-one, period-two, intermittent and chaotic oscillations widely observed in 
experiments. Possible routes to chaos are discussed. 
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§1 Introduction 

The formation of drops is an intriguing phenomenon widely observed in everyday life. Although 
scientific researches on this subject date back to the seventeenth century ,0) great progress has been 
achieved only recently, mainly in detailed studies on the behavior of drops near the breakup point. 



Breakup of a drop is a critical phenomenon corresponding to a singularity of a nonlinear partial 
differential equation obeyed by the fluid with free surface. Refined experiments to observe drops 
falling from nozzles for various viscosities provide much knowledge on this critical phenomenon .SI* 
Eggers proposed a scaling theory which is universally applicable for axisymmetric drops with finite 
viscosity .& Numerical solutions of Navier-Stokes equations have been obtained with high precision for 
viscid0i'i) and inviscid@' fluids, and well reproduce the observed shape of drops. 

Another interesting aspect of the drop formation is the long-term behavior of a dripping faucet as 
a chaotic dynamical system. Since Shaw's pioneering work!) it has been confirmed in many experi- 
ments that dripping time intervals exhibit various types of periodic and chaotic oscillations including 
intermittency and hysteresis Drastic change of an attractor is induced by a small variation 

of the flow rate, which is a main control parameter of the system. However, theoretical progress in 
this direction is not yet satisfactory. A simple "mass-on-a-spring" model has been proposed and an 
analog simulation of this model was reported to reproduce, in a qualitative sense, some experimental 
return maps-Hi) A numerical simulation based on a stochastic model was presented recently, in which 
an Ising-like system was used.© Since the latter model ignores kinetic terms, it might be applicable 
only for very small diameter of faucet(< 0.5mm). In contrast, the mass-spring model corresponds 

* E-mail: fuchi@phys.metro-u.ac.jp 
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to relatively large faucets, where the oscillation of the center of mass governs the basic dynamics. 
However, to explain the complex behavior observed in experiments systematically, information on the 
model parameters (for example, the spring constant) is essential, which is, unfortunately, not easily 
obtained. Both theoretical studies do not aim at investigating the realistic shape of drops, focusing 
on the long-term behavior of the dripping time intervals. 

In our views, the reproduction of the shape of a drop is crucial in general understanding of the 
physics of the phenomenon. Indeed, our detailed simulation turns out to be successful in this respect. 
Even when one attempts to analyze the phenomenon using a simplified model, the knowledge of the 
drop shape is necessary for the correct choice of parameters. Our first goal is a long-term simulation 
using an algorithm which can also simulate the shape of drops reliably. After that, on the basis of 
detailed analysis of numerical data thus obtained, we are going to construct an improved mass-spring 
model which we expect reveals essential features of this complex system. 

In this paper, we present a new algorithm for simulating a dripping faucet based on Lagrangian 
description instead of Eulerian one like Navier-Stokes equations. We decompose a drop into many 
parts (at most 300 ~ 400 disks or less) and describe the dynamics in terms of time evolution equations 
obeyed by each part under the influence of gravity, surface tension and viscosity. After we observe 
that our algorithm can well reproduce shape of a drop at various stages of evolution, we proceed 
to a long-term simulation in which the growth and fission of a drop take place many times under 
a constant increase of the mass of the fluid. It will be shown that time series of the dripping time 
intervals reproduce various types of motions including period-one, period-two, intermittent oscillations 
etc. besides chaotic ones, as observed in many experiments. A bifurcation diagram will be presented 
which agrees well with a recent experiment. 

In §2, we derive static shapes which are used as an initial condition for solving dynamical equa- 
tions. In §3, equations of motion in Lagrangian description is presented. Computational results and 
experimental data are compared in §4. In Appendix A, we describe a variational algorithm to examine 
stability of static solutions. Our algorithm of simulation of dynamics is presented in Appendix B. 

§2 Equilibrium States 

We first derive static equilibrium shapes of a pendant drop where the gravitational and the surface 
tension forces balance each other. There is a maximum (i.e., critical) volume V c for the stable state 
when the radius a of the faucet is fixed. We will see below that there exist in general several equilibrium 
shapes when a volume smaller than V c is given. As suggested by Padday and Pittjli^ only one among 
them is stable and realized, which we will show numerically. 
The force balance equation of the drop is 



where z is the vertical coordinate (the positive direction is defined to be downward), P the pressure. 
The density p is assumed to be constant throughout the fluid. The pressure (difference between the 
inside and the outside of the interface) is expressed in terms of the principal radii of curvature of the 
drop surface, R\, R2 as 



P = 



pgz, 



(2.1) 
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where T is the surface tension. Choosing the length, mass, and pressure units as 



\I V /P9 , 



to 



pll 



(2.3) 



( Iq = 0.27cm, too = 0.020g and Po =270dyn/cm 2 for water at 20 °C,) we can set g 
For an axisymmetric drop, the radii of curvatures are given by 
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where the variables s, r and 9 are defined in Fig. |]. The shape of the drop is thus obtained by 
integrating a set of ODE's: 



dr 

dl 

dz 



sm( 
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d9 

ds r 
0, z(0) = P h , 9(0) 



(2.6) 



tt/2 at s = 0, where P^ is the pressure at the 



with an 'initial' condition: r(0) 
bottom of the drop. 

For a given value of P^, the solution of eq. ( [2.6D is unique. However, we need a boundary condition at 
the top of the drop to determine the upper limit of integration. Here it should be noted that eqs. ( |2.1| ) 
and ( |2.2j ) describe only the force balance, namely, they correspond to an equilibrium condition but not 
to any stability condition. We therefore have to select stable shapes out of the solutions of eq. (|2.6|) . 
It is generally a difficult mathematical problem to examine whether an equilibrium state is stable or 
not when the system energy has not a lower limit, just like the present hanging drop. To the best of 
our knowledge, the existence of a stable shape and its uniqueness have not been proved. Padday and 
Pitt investigated stable shapes of drops under various conditions, assuming implicitly the uniqueness 
of the stable solution. In this paper, we numerically confirm that when there is one or more solutions 
of eq. ( |2.6D , all except one are unstable. The rest is stable under any infinitesimal deformation, which 
we cannot prove but we conjecture. 



2. 1 Drop hanging from a faucet 

If we assume that the inner surface of the faucet is fully wetted and the outer surface is not wetted at 
all, the boundary condition is volume-radius limited@ That is, the shape of the drop is obtained by 
integrating the above equations up to r = a, where o is the inner diameter of the faucet. Solution r of 
eq. (|2.6| ) oscillates as s (and hence z) varies, which leads to, in general, many equilibrium shapes for 
fixed boundary condition. We can show, however, that at most one among them (the shortest shape) 
is stable. 

In the following discussion, we fix at a = 0.5. Figure [2] shows various profiles of drop obtained from 
eq. (|2.6|), The volume of drop, V and its potential energy, U are plotted against the bottom pressure, 
Pb in Figs. ||(a) and[3|(b), respectively. The energy U is the sum of the gravitational energy U g and 
the energy of surface tension Ur- U = U g + Ur- These quantities are calculated from solutions of 
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eq. (2.6). Figure ||(a) indicates that, as a function of V, Pb is generally multi- valued (when V>0.6). 



For example, two possible profiles with the same volume V = 0.6 (having Pb = 3.86 and 4.80,) are 
presented in the third column of Fig. ||[ Also, one can see from Fig. |3|(a) that there are many shapes 
with the same volume, say V = 1.5. Each column of Fig. || presents several equilibrium profiles 
corresponding to a fixed volume. The drop on the top of each column has the lowest bottom pressure 
Pb, an d the second (third) one has the second (third) lowest Pb. These profiles show that when the 
volume is fixed, the lower Pb corresponds to the shorter drop length (= the pressure difference between 
the top and the bottom of drop) as can easily be understood from eq. ( |2.2[ ). 

Observing Figs. ||(a) and ||(b) together, one can see that when the volume is fixed and if there are 
more than one equilibrium shapes (which takes place for V >0.6), the shape with smaller Pb, namely, 
a shorter shape, has a lower energy. This means that a decrease of the surface energy due to lowering 
Pb overwhelms an increase of the gravitational energy induced by lifting up the center of mass. The 
relation such that the lower Pb corresponds to the lower U generally holds in volume-radius limited 
systems but not for volume-angle limited ones like a drop hanging from an infinite horizontal plane. 
(The volume-angle limited systems will be mentioned later on.) Now one might expect that the stable 
shape has the smallest equilibrium energy. But that is not self-evident. In fact, it turns out that the 
statement is true for the volume-radius limited system, but not for the volume-angle limited one. In 
order for any equilibrium shape to be realized, U must increase under any infinitesimal deformation 
with the constraint such that both V and a are fixed. 

In Appendix A, we derive a sufficient condition for stability. Stability of shapes that do not satisfy 
this sufficient condition are numerically examined and it was found that any equilibrium shape except 
the one having the lowest Pb (and hence the shortest one) is unstable under some axisymmetric 
deformation. That was examined for various values of the top radius a. On the other hand, it is 
difficult to determine whether or not the shortest shape, which has the lowest energy, is really stable 
under any infinitesimal deformation. We have studied several types of axisymmetric deformations and 
conjecture that the shortest shape is stable at least under any axisymmetric deformation, so that the 
profile on the top of each column is realized when the volume is given as indicated. The drop on 
the right end in the first row, having a neck close to the tip of the faucet, is also stable. Its volume 
V = 2.39 is almost equal to the critical volume V c indicated in Fig. ||(a). After all, as the volume 
increases, the profile of static drop changes from left to right in the first row of Fig. [2|. 

Interestingly, there are equilibrium shapes with infinite numbers of necks (for V ~ 1.5 when a = 0.5), 
although such shapes are not realized because they are unstable under infinitesimal deformations. 

2.2 Drop hanging from an infinite horizontal plane 

Drops seen for example on the ceiling in the bath-room are volume-angle limited systems. They can 
also be obtained by integrating eq. ( ^ ) up to 6 = ir/2 if the plane is fully wetted. Several equilibrium 
shapes with fixed volume are presented in Fig. || Plots of V vs. Pb and U vs Pb are given in Figs. ||(a) 
and ||(b), respectively. In each column of Fig. B profiles are arranged in increasing order of Pb (and 
hence increasing order of drop length). From a variational analysis, all profiles except the first row 
are found to be unstable under some axisymmetric deformation. As the volume increases, the realized 
static shape of drop changes from left to right in the top row. 

It should be noted that if there are more than one equilibrium shapes, the stable one is the shortest 
but its energy is the second lowest as can be seen from Figs. |](a) and ||(b). This is not surprising 
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because the equilibrium state with the lowest energy is, as pointed out already, not necessarily stable 
for systems whose energy has not a lower bound. In the present case, each shape in the second top 
row of Fig. m has a lower energy than in the top row, which implies that a decrease of the gravitational 
energy overwhelms an increase of the surface energy. But the shape cannot be stable under some 
axisymmetric deformation. 

The critical volume is V c ~ 18.98(= 0.37cm 3 for water at 20° C). A drop larger than this volume 
cannot be suspended from the ceiling. 



§3 Lagrangian Description of Fluid Motion 

We take the following assumptions: 

(1) Incompressible fluid. 

(2) Axisymmetry. 

(3) Horizontal component of the fluid velocity can be neglected in comparison with vertical one. 

(4) Vertical component of the velocity depends only on the vertical coordinate. 

(5) No exchange between upper and lower parts of fluid. 

Assumption (||) is derived from assumptions ( ph a nd (Q). These assumptions correspond to the shallow- 
water theory applied to axisymmetric systems^ and have been used widely. 

Let us denote the volume of the fluid below the vertical coordinate z (the positive direction is defined 
to be downward) by 

t(z,t)= ^r(C,t) 2 dC, (3.1) 

J z 

where z^t) is the vertical coordinate of the bottom at time t and r(z,t) is the radius of the fluid at 
the coordinate z. Then, from assumptions (||), £ can be used as a Lagrangian variable. (See Fig. g.) 
The kinetic energy is thus given as 



#kin - § , 
2 Jo 



/ v 2 dti, (3.2) 
Jo 

dz(£,t) 



v = v z 



at 

Here, the coordinate of the top of the drop (= the end of the faucet) is defined as z = 0, and 
£o(i) = £(0,i) is the total volume of the drop. The gravitational energy is 

U g = -pg z(£,t)d£. (3.3) 
Jo 

The surface energy is expressed as 

C/ r = ry 2irr(z,t)\/l + {dr(z,t)/dz} 2 dz, (3.4) 

which can be rewritten as 



£o(t) j A t i ( z ")2 

" ' (7 



?7r = r/ ^ttz' + Vt^, (3.5) 
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where 

,_dz{Z,t) „_d 2 z(^t) 



Z = 7TZ , Z 



Equations (|3.2|), (|3.3|) and ( |3.5| ) yield Lagrangian of the system as 



C = E kin -U g -U T . (3.6) 

The effect of the viscosity is expressed by a dissipation function, E k - m (namely, the time derivative of 
the Kinetic energy of fluid) asS 

2 V J {<).,; ih-J (37) 

= -3,f (t) /^M} 2 de. 



I 9z 



In the above, we have used the relation 



dv x dv v dv z 

— - H H = . 

dx dy dz 

so that 

OV x OVy 1 dv z 

dx dy 2 dz ' 
because of assumptions ([]]) and (|2|). Equation ( |3.7| ) reads 



-Ekin = — 377 / , , -t , 12 d£ . (3. 
Jo {dz(Z,t)/d£} 2 



We discretize the integral forms ( |3.2| ),(p7^),( |3.5[ ) and ( |3,8[ ). Let the drop be sliced into M disks by 
(M— 1) horizontal planes = Z\, Z2, • • • £m-1 a s Fig. |7[ The volumes of these disks, A£j, are expressed 
as 

A& = P vrr(C,t) 2 dC; i = l,2,...Af; (3.9) 



z A/ = z b . (3.10) 

These volumes are conserved because of the incompressibility. (Variables at the top and the bottom 
of the fluid are defined suitably by taking account of the boundary condition. See below.) The mass 
of the disks are 

Arrij = pA£j , 

which reads 

Arrij = A£j (3.11) 

by taking units such that p = 1. Let the coordinate of the center of mass of the j-th disk be Zj. Then 
the kinetic and potential energies are given by 

^ M 

E kin ~-^Am i (l J ) 2 , (3.12) 
j=i 

M 

U 9 - -d^^mjZj. (3.13) 
3=1 
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The dissipation function is 



where 



^^-^E ^ri H ( 3 - 14 ) 

• =1 \ z j z j~l) 



d 

V i = Z i = W 

The most important part of the algorithm is how to approximate the surface energy because the 
rigorous expression of the surface tension force includes higher-order derivatives in space coordinates. 
We approximate the surface energy U-p = TS (S is the surface area of the fluid) by patching many 
parts of conical surfaces: By defining the average radii of disks as 



Am,- 

— ' j = 2,3,---M; 



3 y n(zj - Zj-x) ' 

the surface area in the interval [(zj-x + z j)/^i i z j + z j+i)/'A is approximated as 



Sj = Tr( rj + r j+ x)\l -Az-j+x + Zj-x? + (rj - r j+1 f 



Sj( z j— l) z jj z j+i) ■ 



Then the surface energy becomes 



U r = T^2Sj(zj-x,Zj,Zj+i) 

3=1 



Further, using an approximation 



Zj + Zj-x 



Zj + Zj-x 



(3.15) 



(3.16) 



(3.17) 



J 2 ' J 2 

and substituting eq. (pTT7|) into ( g3g ) and ( gig ), we obtain a Lagrangian of the discretized system as 

£( zx ,---z m ,zi,---zm) 

= ^kinC^l, ' ' ■ z m) -U g (zx,--- z m) ~ U V {zx, ■ ■ ■ Z M ) ■ 



Equations of motion in the Lagrangian description are thus obtained from 

d dC dC 1 dEkin 



dt dzj dzj 2 dij 



; j = i,2,---m. 



In the present simulation, we used a further approximation 



Zj ~ Zj, 



(3.18) 



(3.19) 



(3.20) 



in place of eq. ( 3.17 ). 

The top boundary is treated as follows: We mark a part of the fluid just inside the faucet and let 
its coordinate be zq(< ). The volume in the interval [zq, zx] is 



21 



ACx= irr((,tyd( = 7ra 2 \z \+ Trr((,tyd(, 



Zl 



(3.21) 
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which is constant as time passes. We fix the velocity of the fluid at the end of the faucet so that 

zo = v ; 

namely, the volume of the fluid hanging from the faucet increases steadily with the constant flow rare 
Q = ira 2 vo. The volume of the first disk i.e., the volume in the interval [0, zi], A£ 1; increases as 

A£l = f * vrr(C,t) 2 dC = - Tra 2 \z \ • 
Jo 

The coordinate z\{> 0) increases according to the above relation. When A^ reaches a certain constant 
value, we redefine the numbers as 

i + l<-j; j = l,2,---M; (3.22) 

and reset zq and z\ to be 

zq < A^i/ira 2 , z\ <— . 

In this way the number of disks M increases by one. M also increases when the relative thickness of 
any disk exceeds a certain limit and is divided in two. A more detailed description of the algorithm 
will be given in Appendix B. 

§4 Results and Discussions 

We choose to = (T/pg 3 ) 1 ^ as the time unit. For water at 20 °C, to = 0.017 s. The unit of viscosity 
was chosen as 770 = (/?r 3 / g) 1//4 . Then viscosity is n = 0.002 for water at 20 °C. (Other units have been 
defined as eq. Q2.3| ).) Setting values for the faucet radius a and the velocity vq, we simulate the time 
evolution of the dripping faucet. 

4-1 Comparison of shapes with experiments 

Peregrine et al. observed details of the shapes of drops falling from a capillary tube.i* We first see 
that the present algorithm can well reproduce their experimental data. The faucet radius was chosen 
to be a = 0.952 corresponding 5.2 mm in diameter in their experiment. The flow rate of the water is 
not mentioned explicitly except their expression 'as slow as possible'. We chose as vo = 0.01(= 0.16 
cm/s). This is the velocity at the top of the drop (i.e., the tip of the faucet) employed as a boundary 
conditions. The breakup parameter is e = 10 -4 , which means when the minimum cross section of 
the necking region divided by the cross section of the faucet, S m i n = (r m j n /a) 2 , reaches this value we 
separate the drop from the remaining part and continue the simulation for the residue. 

In Fig. ^, we present profiles of drop falling from the faucet at various stages of evolution. The 
initial shape is taken from a static stable state of the volume V- m \t = 4.77, which has been obtained by 
integrating eq. ( |2.6j ) with = 2.6. The time t = 12.57 is the breakup moment, namely the moment 
at which S m i n reaches e. The time t = 12.67 is the second breakup moment. In Fig. ||, the profiles 
near the breakup moment are superimposed with the experimental photographs. 

We see that the present simulation reproduces the observed shapes excellently. Especially, our shape 
of the satellite (the secondary small droplet) is closer to the observed one than that obtained from a 
different algorithm in which the viscosity is ignored^ (The shape at the first breakup moment has 
been well reproduced by Eggers.0)) If we remember several approximations made in our algorithm 
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(particularly, eq. ( 3.20| ) and the procedure (12) in Appendix B), the agreement between simulation 



and experiment is amazingly good. 

We can estimate an error of the breakup moment based on the following scaling relationship 

r(t) oc (t c -t) 2/3 for inviscid fluid, (4.1) 
r(t) oc (tc-t) 1 for viscid fluid, (4.2) 



where t c is the critical point. The scaling region where eq. (4.2) holds is extremely small for water 



due to small viscosity (typical scales are l v = rj 2 / pT ~ 1.4 x 10 _6 cm, t v = rj 3 / pT 2 ~ 1.9 x 10~ 10 s). On 
the other hand, eq.( [4.1[ ) holds only approximately because of a crossover effect.!^ Nevertheless, these 
relations are useful enough to estimate the upper and lower limits of t c . Using eq. (4.1) and eq. (4.2) 



to extrapolate the plot of r(t), we obtained , for example, 12.58 < t c < 12.66 corresponding to the 
first breakup moment t = 12.57. 

4-2 General feature of drops falling from a faucet 

As has already been recognized experimentally^ and theoretically liquid drops do not have up- 
down symmetry like Fig. |l^(a) near the critical point at which a drop separates. That is not due to 
the gravity. Instead, the symmetry spontaneously breaks down as sketched in Fig. [i~0|(b) or (c) even 
without the gravitational effect. (The shape (c) corresponds to breakup of a satellite.) As we have 
seen, the marked asymmetry is observed in the present simulation as well. It should be noted that the 
asymmetry is a dynamical property of the fluid with free surface, therefore it cannot be predictable 
from the equilibrium shapes as in Fig. |2| and also not reproduced in the simulation neglecting kinetic 
terms.lll' 

Theoretically, from scaling properties derived by Eggers,B0) the asymmetric shapes near the critical 
point are expected to be universal irrespective of the viscosity, the faucet radius and the flow rate 
except the inviscid case. Practically however, it is also well known that the global shapes are strongly 
dependent on these conditions because the scaling law, even if it is exact, works only in the small 
region of time and space close to the critical point. In fact, the liquid bridge that connects the conical 
part (just below the faucet) and the spherical part shows a rough tendency to become long and thin 
for (1) large n, (2) large a and (3) small vq. These features are confirmed also in our simulation. For 



example, we present a simulation for large n. Figure [llj shows a temporal change of a drop of glycerol 
corresponding to an experiment by Shi et al.& Similar profiles have been obtained by solving the 
Navier-Stokes equations.!! Both results well reproduce a long bridge observed in the experiment. 

Strictly speaking, the length of the liquid bridge does not change monotonically with these parame- 
ters because of spatial and temporal oscillations of the liquid surface. The necking region is delicately 
affected by the interplay of these oscillations and liquid influx. Concerning the third condition (i.e., 
^o-dependence), the shape tends to a final one in the limit of vq — > which should coincides with 
that obtained by setting vq = and starting from an initial state with volume just above the critical 
value V c . For example, in case of parameter values (rj = 0.002, a = 0.952, v = 0.01, e = 10" 4 ), 
with the initial volume VJnit = 4.77, the first breakup occurred at the drop size = 3.85 and the 
residue V r = 1.28 (Fig. |i~2|(a), which is the same profile at t = 12.57 in Fig. |8|). These results are 
to be compared with Vj = 3.72 and V T = 1.23, which were obtained by a simulation with no flow: 
vq = 0, the same values for other parameters and the initial volume VJ n it = 4.95 slightly larger than 
the critical value (Fig. |i~2](b)). Two profiles look almost the same. 
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4-3 Formation of secondary drop 



In the experiments focusing on long-term behavior of dripping faucets, the flow rate is usually chosen 
as a control parameter. As mentioned above, the shape of drop is strongly affected by the velocity 
vq or, equivalently, the flow rate. In Figs. [13(a) and |l3|(b), we present how the shape of drop at the 
second breakup moment depends on vq. Parameters are chosen as (a = 1.0, n = 0.002, e = 10" 3 ) 
commonly in (a) and (b); vq = 0.003 in (a) and vq = 0.3 in (b). The initial shape was commonly taken 
to be a static stable one. In both (a) and (b), two profiles describe the first and the second breakup 
moments (i.e., 5 m j n = e). 

For smaller v$, after the spherical drop detached from the bottom of the neck, a satellite droplet 
forms from the slender neck itself and separates from the conical part just below the faucet. The 
satellite is thus much smaller than the main drop. The volume ratio is 0.7 % in the present example. 
For larger vo, the volume increases so rapidly that the slender liquid bridge disappears, the spherical 
part and the cone jointing directly. In the latter case, the secondary drop is again spherical and 
relatively large, its size being 44 % of the main drop. Results for the secondary-drop formation similar 
to Fig. 13 have been obtained from the simulation of inviscid fluid.U 

Now a question is whether such successive processes, breakup of a main large drop followed by a 
smaller secondary drop can take place regularly or not. Details of shapes after the second breakup 
moment has not been reported so far. Figure 14 represents a continuation of Fig. 13(b), namely, a 
series of breakup moments for the same parameter values: a = 1.0, vq = 0.3, where the first profile 
represents the third breakup moment. There appears sometimes a tiny droplet (smaller than ~ 0.01 
% of the largest drops) splitting from the tip of the liquid cone just after the main drip, that is not 
presented here. Computational error may have led to such small droplets but the possibility that 
extremely small droplets really appear cannot be excluded. A larger drop and a smaller drop appear 
alternatively in Fig. 14, which looks almost periodic. Note, however, that large (small) drops are 
slightly different from each other in size and shape. These fluctuations are surely intrinsic, not due 
to computational error. A minute difference at each breakup moment can be an origin of not only 
fluctuations but also an irregular motion under certain conditions. Here we see that the dripping 
faucet is really a complex system. 



4-4 Long-term behavior 

We present here results for o = 0.916 corresponding to the faucet of 5mm in diameter used in a recent 
experiment by Katsuyama and Nagata.0) For the long-term simulation, a larger value for the breakup 
parameter e = 4 x 10 -3 was used to save the computational time. The velocity vq is chosen to be 
relatively small (< 0.16), so that the secondary drop formation is like Fig. |l3|(a) instead of Fig. |i~3|(b). 
In other words, pinching off of a satellite (or occasionally two successive satellites) smaller than 1 % 
of the largest drops always occurs just after breakup of a main drop. In the following analysis, we 
ignore these satellites. Then, the drop sizes distribute mostly in the range larger than 10 % of the 
largest drops. 

We let T n denote the n-th dripping interval, i.e., the time difference between the (n + l)-th and n-th 
drips. Neglect of satellites above mentioned implies that intervals T n which are smaller than ~ 0.2 
are omitted, while T n for main drops are larger than ~ 10. Figure |l5|(a) represents plotting of T n vs 
n for various values of vq. Corresponding return maps, i.e., plots of T n+ \ vs T n and power spectra are 
given in Figs. [l5| (b), and (c), respectively. 
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The time series data {T n } fluctuate considerably, but the power spectra can suggest types of bifur- 
cation. As the control parameter vo is varied, period one (PI) motion at vq = 0.083 period-doubles 
backward to period two (P2) motion. That can be confirmed from the power spectrum for uo=0.0825 
which exhibits a peak at frequency / = 1/2, indicating P2 oscillation. When vq = 0.082, intermittent 
period three (P3) motion is observed, where the spectrum has a sharp peak at / = 1/3. Another type 
of P2 motion is observed at vq = 0.0815. We found that the bifurcation to this type of P2 motion 
is not period doubling from any PI motion. Rather, forward and backward period doubling cascades 
to chaos, starting from this type of P2 motion, are expected. In fact, the spectra for vq = 0.0813 
and 0.0809 exhibit a peak at / = 1/4 besides a sharp peak at / = 1/2, which indicates that the P2 
motion at vq = 0.0815 period-doubles backward. The backward period doubling is widely observed 



also experimentally. Attractors in Fig. 15(b) closely resemble experimental ones. 



Figure |16| is a bifurcation diagram in which the parameter range of Fig. |15| is indicated. As sug- 



gested by Fig. 16, oscillations like those shown in Fig. [D] are quite typical in the present system, and 
repeatedly appear at different values of vq. Similar types of oscillations are observed not only in the 
experiment by Katsuyama and Nagata for the same faucet size,@* but also in other experiments with 
smaller faucets. 

In the present simulation, at vq = 0.074 and vq = 0.083 for instance, the T n value is uniqe. In 
contrast, between these v$ values (for 0.074 < vq < 0.083 for instance), the T n value is distributed 
over a finit range. This is seen on Fig. fill in the form of a bloc, which repeats as vq increases. We may 



call this pattern a unit structure. Looking at Fig. 16 with reference to T n , one sees each unit structure 



occurs periodically, namely, at the same interval in T n . Qualitative agreement of our bifurcation 
diagram with the experiment by Katsuyama and Nagata is satisfactory in a wide range of the control 
parameter vq, or equivalently the flow rate. In their bifurcation diagram, a unit structure similar to 
ours also appears repeatedly as the flow rate is varied. 

4-5 Further study 

We have seen that the present simulation reproduces various aspects of the dripping faucet system 
which are observed experimentally. The next question is why the dripping faucet behaves like this. 
We will answer to this question in the forthcoming paper. Here we give only an outline of it: 

We construct an improved mass-spring model based on a detailed analysis of our numerical simula- 
tion. The model reveals the basic mechanism of the complex behavior of the dripping faucet system. 
A key of the new model is that the mass dependence of the spring-constant is taken into account. A 
similar global feature observed in the experiment and the simulation, i.e., the repeating unit structure 
in the bifurcation diagram, is reproduced also in the improved mass-spring model. This simplified 
model shows that each unit in the bifurcation diagram includes ordinary period doubling cascade to 
chaos, forward and backward period doubling cascade starting from P2 motion, intermittency and 
hysteresis. It is also clarified that the unit structure can be explained in terms of oscillations of the 
center of mass of the fluid, in other words, each unit is characterized by an integer which relates to 
the frequency of the oscillations during the time between two successive breakup moments. 
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Appendix A: Stability Analysis 

We consider an axisymmetric deformation such that the part between the interval [z,z + dz] is 
mapped onto the interval [Z, Z + dZ], where 

dZ= [l + e(z)]dz 

with the constraint 

e(0) = e(*b) = 0. (A-l) 

The coordinates at the top and the bottom of the drop are defined as z = and z^, respectively, which 
are mapped onto Z = and Z^ by the deformation. 

Because the volume of each part is conserved before and after the deformation, the radius r(z) is 
transformed as 

R{z) -^frWy 

The increment of the surface energy Uy caused by the deformation is expressed as 



5U r = 2vrr J " Ryfl + (R') 2 dZ 
- 2irT J h r^l + (r') 2 dz 



= 5U r ,i + SU r ,2 , 

where 8Ut,i = 0(e) and 5Ur,2 = 0(e 2 ) are the first and the second order quantities of the small 
deformation e. The increment of the gravitational energy U g is 

f Z b „ f«b „ 

5U g = —pg l ttR ZdZ + pg / nr zdz = 5U g i; 
Jo Jo 

5U g includes only linear terms of e, namely the second order term vanishes: 5U gt 2 = 0. It is a little 
tedious but not difficult to derive the force-balance equation fl2.ip and ( [2. 2D from the equilibrium 
condition 

SU T ,i + 6U gA = 0. 



The stability condition is 
where 

SU r ,2 

Hz) 



SUr 2 > 0, 



r (z)e 2 + iP(z)(e') 2 ]dz, (A-2) 



o 



4 

(T+ (r') 2 ) 5/2 
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+ (r'f + (r') A -rr" (j -2(r') 2 ') 



(A-3) 



«"> - V^W (A ' 4) 

The second term in the integrand is positive definite. Therefore the condition 

<f>(z) > for < z < z h (A- 5) 
is sufficient for the stability of equilibrium shape under any axisymmetric deformation e(z). 



When the drop has a neck, the condition ( A-5 ) is not satisfied because (j>(z) is negative at the neck 
(r' = 0, r" > 0). We have found that all shapes with volume V > 1.3 in Fig. |2| and V > 17.0 in Fig. || 
violate this sufficient condition. (Some of them have no neck but do not satisfy ( |A-5| ).) However, such 
shapes are not necessarily unstable because the the homogeneous deformation e 1 = is impossible 
under the constraint of eq. ( |A-1[ ), and hence the second term of the integrant in eq. ( |A-3| ) always 
works in the direction of stabilization. 

Therefore we first calculate (j)(z) and if 4>{z) < in an interval < z\ < z < z\ + d, the following 
infinitesimal deformation is considered: 

e{z) = 2A {Z ~ Z2) 



e(z) = 2A 



(ad) 2 (3 

for Z2 < z < f3ad + Z2, (A-6) 
(ad + Z2 — z) 



(ad) 2 (l- (3) 

for j3ad + Z2 < z < ad + Z2, (A-7) 
e(z) = otherwise, (A-8) 



which leads to a stretch of the original shape by the length 

fZ2+ad 

edz = A . 



The parameters Z2 (satisfying < Z2 < Z2 + ad < z\,), a (> 0), and (3 (> 0) characterizing the 
deformation are varied and we examined if SUy^/^ 2 becomes negative. In this way, for fixed values 
of V and a, all the equilibrium shapes except the shortest one, which has the lowest energy, are found 
to be unstable for the volume-radius limited systems. Also for the volume-angle limited systems, all 
the equilibrium shapes except the shortest one, which has the second lowest energy, are found to be 
unstable. All the shapes (including the largest one) on the top of each column in Figs. ^ and || are 
found to be stable at least under axisymmetric deformations. 

Appendix B: Algorithm 

We employed a fourth-order Runge-Kutta method with adaptive stepsize control. When the 
radius a of the faucet and the velocity t>o of liquid at the tip of the faucet are given, the simulation is 
done based on the following algorithm. 

(1) Initial shape 

We choose one of the stable equilibrium shapes obtained from the method described in r2. 
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(2) Decomposition of a drop 

The drop of (jlj) is sliced with horizontal planes so that the drop is decomposed into many disks. 
The thickness of disks are such that the length s of each disk measured along the edge of the 
section shown in Fig. |7] is equal. Therefore the disks near the bottom are relatively thin(see 
Fig. fj). The total number of the disks is denoted as M. 

(3) Coordinates and volumes of disks 

Vertical coordinates of the slice plane are denoted as zi, z%, ■ ■ ■ zm-i- They satisfy the relation 

< z x < z 2 < ■■■ < z M = Zb, (BT) 

where 

Zl (t = 0) = 0. 

Further, we set zq to be a negative constant value: zq = Zo(t = 0) < 0. The volume of the disk 
in the interval [ denoted as A£,-(j = 1, 2, • • • M) and the volume in the interval [0, z\] 

as A^ x . These volumes are calculated precisely from the static solution employed as an initial 
condition. 

(4) Average radius of disks 

We define the average radii of disks as follows: 



A6 



ir(zj - Zj-i) 



j = 2,3,---M; 



n = \ (zi>0) 



TTZl 

(5) Initial velocities of disks 

We set as Vj(t = 0) = vq (j = 1, 2, • • • vm)- 

(6) Tentative stepsize 

At t = 0, we tentatively set as At = At. 

(7) Adaptive stepsize 

From the time-evolution equations 



(B-2) 



fj({zih {n(zi,zi-i)}, {vi}) 
i = l,2,...M; 



(B-3) 



we estimate the error A for one step of fourth-order Runge-Kutta algorithm. Then A is compared 
with our desired accuracy An and an adjusted stepsize At is decided. This procedure is applied 
for all errors of variables { 2i} , („,} and the smallest value of At is take,,B 

(8) Time evolution 

Equation ( |B-3j ) is integrated by one time step i.e., from t to t + At using the fourth-order 
Runge-Kutta method. 
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(9) Renewal of the stepsize 

From (|8|), we estimate the error A in a similar way to (|7|) and renew the adapted stepsize At for 
the next step. 

(10) Check of overtaking 

Because of assumption (||) in §4, any disk must not overtake the neighboring one, in other words, 
the relation ( |B-1| ) should be maintained all times. So if the condition of eq. ( |B1| ) is violated, 
we choose a smaller value for At and redo the same procedure (||). This does not occur if the 
accuracy Ao has been chosen small enough. 

(11) Increment of the number of disks 

When the volume of the first disk, A^ becomes larger than a certain constant value, we redefine 
variables as 

Zj+i <- Zj\ j = l,2,---M; 
Vj+i <- Vj\ j = 0,1, ■■■M; 
zi <- 0, 

and reset zq to be the initial value Zo(t = 0). In this way, M increases by one. 

(12) Division of a disk 

If the ratio of the width to the radius, {zj — Zj_i)/rj is over a certain limit for any disk, the 
disk is divided into two parts so as to conserve both the total volume and the momentum. In 
this way, M is increased by one. It should be noted that the way of division is not unique. On 
the other hand, it turns out that there is no solution of division under additional constraint 
such that the energy is rigorously conserved. We divided disks in such a way that the spatial 
derivatives of velocity is conserved besides the volume and the momentum. 

(13) Coalescence of disks 

When the radii of some neighboring disks become larger than a certain value, the two disks are 
combined into one. Then M decreases by one. 

(14) Renewal of disk radii 

From the renewed values of {zi} and {A£j}, the values of {r,} are replaced using eq. (|B-2| ). 

(15) Breakup of the drop 

When the shape of the drop has a neck (or necks), the minimum radius r m [ n in the necking 
region is compared with the faucet radius. If its relative cross section S m m = (r m /a) 2 is less 
than a critical value, then the part below the neck is separated. After that, the number of disks 
becomes M <— (m — 1) where in is the disk number corresponding to r m — ^min* 

(16) Return to (§) 

We go back to the procedure (|8|) and continue the iterations. 
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Fig. 1. Definition of variables. 
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Fig. 2. Drop hanging from the faucet, a = 0.5. Various equilibrium shapes for a fixed volume of the 
drop. Among them, only one (the top of each column) is stable. As the volume is increased, the shape 
changes as indicated by arrows. The maximum volume is V = 2.39. 
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Fig. 4. Same as in Fig. § but for a drop hanging from the ceiling. 
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/ \- z = z&t) 

V * J 

' — •-• z-z h {L) 

Fig. 6. Definition of variables. 
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1-- z = 


A^ 2 \ 



z = z (t) 
z = z x {t) 
z = z 2 (t) 



z=z M _At) 



z = z M {t) 



Fig. 7. Definition of variables. 




Fig. 8. Temporal change of the shape, a = 0.952, v$ = 0.01, rj = 0.002, e = 10 4 . Initial condition: 
Pb = 2.6 corresponding to the initial volume Vj n i t = 4.77. 
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Fig. 9. A comparison of experimental and calculated drop shapes close to breakup moments. The 
photographs of the drops are from Peregrine et al. (J. Fluid Mech.,212 (1990) 25. Reprinted with 
permission of Cambridge University Press.) a = 0.952, v = 0.01, 77 = 0.002, e = 10~ 4 . (a) The 
moment at which S m i n = 5e. (b) The moment at which S m i n = e. 



(a) (b) (c) 




Fig. 10. The shape at the critical point is like (b) or (c) instead of (a), (c) is for pinching of a satellite. 
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t = 12.8 



t = 13.8 



t = 14.5 



Fig. 11. Time evolution of a drop of glycerol. The length, time and viscosity units are = 0.266 cm, 
to = 0.0125 s and rjo =4.26 g/cm s, respectively, corresponding to glycerol at 20 °C. a = 0.332 (1.5 
mm diameter), r\ = 3.50, vq = 0.05. 





Fig. 12. Profile at the critical point in the limit of v -> 0. a = 0.952, 77 = 0.002, e = 10" 4 . (a) 
v = 0.01. (b) v = 0. 
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Fig. 13. Formation of the secondary drop, a = 1.0, 7] = 0.002, e = 10~ 3 . Initial condition = 2.6 
corresponding to Vinit = 5.21. Profiles are at the breakup moment, at which .Smm = e. The time of 
the breakup moment is written in each frame, (a) vq = 0.003. (b) vq = 0.3. 
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Fig. 15. (a) Plot of dripping time interval T n vs n for various values of vq. Value of t>o is written in 
each frame, (b) Return map: plot of T n+ i vs T n . (c) Semi-log plot of power spectrum of (a) calculated 
from 2 8 data points. 




Fig. 16. Bifurcation diagram of the dripping interval T n vs vq. 
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